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A correct description of electronic exchange and correlation effects for molecules in contact with extended 
(metal) surfaces is a challenging task for first-principles modeling. In this work we demonstrate the impor¬ 
tance of collective van der Waals dispersion effects beyond the pairwise approximation for organic-inorganic 
systems on the example of atoms, molecules, and nanostructures adsorbed on metals. We use the recently 
developed many-body dispersion (MBD) approach in the context of density-functional theory [Phys. Rev. 
Lett. 108 , 236402 (2012); J. Chem. Phys. 140 , 18A508 (2014)] and assess its ability to correctly describe the 
binding of adsorbates on metal surfaces. We briefly review the MBD method and highlight its similarities to 
quantum-chemical approaches to electron correlation in a quasiparticle picture. In particular, we study the 
binding properties of xenon, 3,4,9,10-perylene-tetracarboxylic acid (PTCDA), and a graphene sheet adsorbed 
on the Ag(lll) surface. Accounting for MBD effects we are able to describe changes in the anisotropic polar¬ 
izability tensor, improve the description of adsorbate vibrations, and correctly capture the adsorbate-surface 
interaction screening. Comparison to other methods and experiment reveals that inclusion of MBD effects 
improves adsorption energies and geometries, by reducing the overbinding typically found in pairwise additive 
dispersion-correction approaches. 


I. INTRODUCTION 

Atoms, molecules, and extended nanostructures in¬ 
teracting with solid surfaces are omnipresent in mod¬ 
ern materials^. Obtaining a detailed understanding of 
the physical and chemical properties of metal-adsorbed 
molecules is of paramount importance in a wide variety 
of fields. However, the description of such complex sys¬ 
tems poses a significant challenge to first-principles mod¬ 
eling by featuring all limiting cases of chemical bond¬ 
ing ranging from the delocalized nearly-free electrons in 
the metal to directional covalent bonds in the adsor¬ 
bate. In addition to these limiting cases, the descrip¬ 
tion of adsorbate-surface bonding and dispersion inter¬ 
actions create a range of additional challenges. Many 
recent works suggest that dispersion interactions play a 
dominant role in many realistic hybrid organic-inorganic 
materials^^®. The interplay between localized and de¬ 
localized states at such an interface can be translated 
into the relevance of both typically distinguished do¬ 
mains of electron correlation: the “static correlation” 
or inherent state-degeneracy that underlies the substrate 
metallic states, as well as the “dynamic correlation” that 
governs the long-range dispersive interactions induced 
by the quantum-mechanical fluctuations within the com¬ 
bined adsorbate-substrate system. Quantitatively or of¬ 
ten even qualitatively correct first-principles treatment 
will have to efficiently account for both of these effects^’®. 
Whereas wavefunction-based approaches to electron cor¬ 
relation yield an excellent description of the latter®^^^, 
density-functional theory (DFT), specifically in its semi¬ 
local approximations based on the homogeneous electron 
gas, is able to describe the metal electronic structure rel¬ 
atively well^^d3_ 

Metallic states can be treated correctly by utilizing 
semi-local approximations to DFT such as the local- 


density approximation (LDA) or the generalized gra¬ 
dient appoximations (GGAs) that satisfy the homoge¬ 
neous electron gas limit^^’^^. On the other hand, long- 
range correlation and dispersion interactions of isolated 
systems is best treated with well-established quantum- 
chemical approaches, such as the coupled cluster (CC) 
method^®. However, long-range correlation effects in 
solids or surfaces such as dispersion interactions between 
adsorbate and substrate, are still out of reach for these 
approaches and completely amiss in semi-local DFT^®’^^. 
The resulting lack of van der Waals (vdW) interac¬ 
tions can lead to failure to find any stable adsorbate 
structures^®d9. Accurate treatment of long-range corre¬ 
lation and their manifestation as dispersion interaction is 
key to correctly describe the bonding of adsorbate-surface 
complexes. While many approaches beyond semi-local 
DFT exist trying to incorporate an improved descrip¬ 
tion of exchange and correlation, either on the basis of 
the adiabatic connection fluctuation-dissipation theorem 
(AGFDT)^®’^^ or many-body perturbation theory^^^^^, it 
may still be desirable to retain the simplicity and compu¬ 
tational efficiency of semi-local functionals, but somehow 
incorporate a physically correct description of long-range 
correlation as it is given by wavefunction approaches. 

Empirical pairwise additive dispersion-correction 
schemes on top of DFT, such as DFT inch dispersion 
(DFT-D) in the variants proposed by Grimme^^’^® can 
be seen as a first step to form such a hybrid approach. 
A simple physical form of the leading-order dispersion 
term is included and the corresponding empirical atomic 
parameters such as dispersion coefficients Cg, atomic po¬ 
larizabilities, and van der Waals radii are precalculated 
and tabulated. A transition to the pure DFT treatment 
is assured by damping the vdW interactions at close 
interatomic distances. A more general and less empirical 
treatment of dispersion interactions is provided by the 
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approaches of Becke and Johnson^^ or Tkatchenko and 
Scheffler (TS)^®, where the dispersion parameters are 
functionals of the electron density calculated from DFT. 
The corresponding parameters therefore adapt to the 
chemical environment of the atom. Due to their success 
and efficiency, pairwise additive treatments of dispersion 
are now standard practice in gas-phase quantum chem¬ 
istry in describing molecular aggregates and also yield 
a qualitatively correct description of metal-adsorbed 
molecule geometry and energetics as has been shown 
in the case of benzene and azobenzene adsorbed on 
coinage-metal surfaces^®’^®. Nevertheless, comparison to 
experiment has shown that a simple pairwise treatment 
of dispersion can induce significant overbinding, hence 
quantitative or “predictive-quality” description cannot 
be achieved in general^®. This is due to the fact that 
many-body correlation effects and Coulomb-screening 
within the substrate play a significant role in these 
systems®®. Accounting for the latter by effective 
renormalization of metal Cq coefficients, as done in the 
vdW®“''^ scheme,®® has improved adsorbate geometries 
significantly®^, albeit at a remaining overbinding in 
adsorption energies. Furthermore, most studies of 
molecules on surfaces have concentrated on the limit of 
low coverage. Obviously, many-body dispersion effects 
will become even more pronounced for dense molecular 
monolayers and multilayers adsorbed on surfaces. 

A different approach to include dispersion interactions 
in DFT is given by the vdW-DF type of density func¬ 
tionals, where an additional non-local correlation con¬ 
tribution is added to a semi-local exchange-correlation 
functional®®’®®. The non-local correlation is described 
by a two-point integral over the electron density and a 
given integration kernel. Although the correlation inte¬ 
gral is of a two-body form, higher-order semi-local con¬ 
tributions can be effectively incorporated in the formu¬ 
lation of the kernel. Recent improvements in computa¬ 
tional efficiency®^ and performance®®’®® have triggered a 
more widespread use of vdW-DF specifically in the con¬ 
text of metal-surface adsorption®®. Problems of incon¬ 
sistent exchange treatment in earlier versions, resulting 
in systematic underbinding, are being addressed in more 
recent versions, such as the vdW-DF-cx functional®^’®®. 

An efficient long-range approach to correlation that 
goes beyond pairwise dispersion has recently been pro¬ 
posed in the form of the many-body dispersion or MBD 
method"*’’®®. In the MBD approach, the problem of cal¬ 
culating the long-range correlation of a set of atoms in a 
molecular arrangement is recast to calculating the corre¬ 
lation energy of a coupled system of quantum harmonic 
oscillators (QHO) interacting via the long-range dipole 
potential. Both, many-body contributions and screen¬ 
ing effects due to mutual polarization of the QHOs, are 
accounted for and a range-separation in terms of the 
Coulomb potential facilitates the connection to the DFT 
functional®. This approach quantitatively describes in¬ 
teraction energies of a wide variety of systems including 
many well established gas-phase testcases"*’^® as well as 


molecular crystals^*^ and supramolecular arrangements'^® 
at minimal overhead compared to the pure DFT calcu¬ 
lations. Furthermore, the foundation of MBD on the 
ACFDT enables a systematic improvement in terms of 
better approximations to the response kernel or using 
well-established quantum-chemical techniques for elec¬ 
tron correlation as applied to coupled oscillators. 

In this work we aim to assess and establish the MBD 
approach for hybrid organic-inorganic systems, specif¬ 
ically organic adsorbates on metal surfaces, a field in 
which MBD contributions to binding play a dominant 
role. After shortly revisiting the method and discussing 
specific issues in the context of metal surfaces, we inves¬ 
tigate the relevance of MBD interactions for an atom, a 
molecule, and an extended nanostructure adsorbed on a 
metal Ag(Ill) surface. We find that in all cases MBD 
interactions are highly important in order to account for 
the correct interaction energy and even more so to cor¬ 
rectly describe response properties such as vibrational 
frequencies or the polarizability of the system. Com¬ 
paring to experiment and other simulation approaches, 
we find significant quantitative improvement by includ¬ 
ing many-body dispersion effects. We conclude the work 
by shortly outlining the remaining challenges to establish 
the many-body dispersion approach as a contender in the 
modeling of molecule-surface systems. 

II. THE MBD METHOD 

In the following we will review the physical foundations 
of the MBD method as it was recently published and ana¬ 
lyzed in detail elsewhere^’®’®®. The quantum-mechanical 
electron correlation energy, which contains the dispersion 
energy of a system, can be calculated from the micro¬ 
scopic density-density response function x(r, r', ito) using 
the ACFDT®’43’44: 

1 r°° 

= d. ( 1 ) 

X / dATr[(xA(r,r',ia;) - xo(r,r',fa;))u(r,r')]. 

Jo 

In this formalism the response function x at a certain 
interaction strength A is calculated self-consistently from 
a non-interacting reference response function xo(r, ioj) 
using a Dyson-like screening equation^®: 

XA = Xo + Xo(Au-k/a;c)XA, (2) 

with V being the Coulomb potential and fxc being 
the exchange-correlation kernel. Approximations to the 
above equations typically vary in the initial guess for 
the non-interacting response function xo and in the ap¬ 
proximations to fxc (e.g. fxc = 0 in the random- 
phase approximation (RPA)). Typically xo is con¬ 
structed using a set of effectively independent particles 
or quasiparticles^®’^^, such as Hartree-Fock (HF) states 
or Kohn-Sham (KS) states as a result of Hartree-Fock or 
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DFT calculations^®’^®. These (quasi)-independent states 
are renormalized with respect to the Coulomb interaction 
via Eq. (2). Approaches of this kind to calculate the cor¬ 
relation energy have recently been established using the 
RPA on top of single-particle states from DFT®®’®®. Even 
in the case of RPA, solving the above equations quickly 
becomes computationally intractable with system sizes 
that are typically needed to model realistic materials. 

The main idea behind the DFT-I-MBD approach is to 
find an alternative efficient way to construct xo solve 
the correlation problem with the help of an intermediate 
set of quasiparticle states. The effects of mutual instanta¬ 
neous polarization and depolarization of (valence) elec¬ 
trons can be recast into a system of effective quantum 
harmonic oscillators (QHOs) or Drude quasiparticles®^. 
The electron density around every atomic nucleus is rep¬ 
resented by such a three-dimensional QHO with an ef¬ 
fective width, mass, and frequency that connect to the 
polarizability and dipolar response of the valence elec¬ 
trons. This so-called coupled fluctuating dipole model 
has proven very successful in describing long-range in¬ 
teraction and polarization®®^®^. In this picture the non¬ 
interacting response xo of a system simplihes to a sim¬ 
ple product of individual localized QHO quasiparticle 
response functions. The corresponding interacting re¬ 
sponse function and correlation energy can then effi¬ 
ciently be calculated within the framework of Eqs.(l) and 
(2) using any approximate quantum chemistry method. 

The challenge of recasting the long-range correlation 
problem to a set of QHOs lies in hnding a connection 
to current electronic structure methods correctly treat¬ 
ing short-range interactions, such as semi-local density- 
functional approximations. Treating the long-range cor¬ 
relation problem with QHO quasiparticles, interactions 
between these should only include terms not yet treated 
in the short-range via DFT. Correspondingly the quasi¬ 
particles need to already effectively contain short-range 
polarization and screening effects that are included in the 
DFT description. One therefore defines the frequency- 
dependent dipole polarizability of every QHO 




apHr)] 

1 -b (w/wp[n(r])2 
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via the static atomic polarizability ap[n(r)] and the 
characteristic excitation frequency a;p[n(r)] of the atom 
it models. These parameters can be extracted from 
the electron density as predicted by DFT using the 
Tkatchenko-Scheffler (TS) scheme®®. Therefore 
includes hybridization effects due to the local environ¬ 
ment as well as short-range exchange-correlation effects. 
Short-range interaction screening and anisotropic polar¬ 
izability changes are furthermore included by renormal¬ 
izing the polarizability with respect to the short-range 
Coulomb interaction of two spherical gaussian charge dis¬ 
tributions associated with each pair of oscillators®. The 
resulting QHOs (described by screened polarizabilities dp 
and screened characteristic frequencies Wp) are quasipar¬ 
ticles that implicitly contain the short-range polarization 


due to the presence of other QHOs. The partition into 
long-range and short-range interactions is made using a 
standard range-separation technique, which, due to the 
arbitrariness of this partition, introduces a single range- 
separation parameter /3, which has to be predetermined 
once for a given exchange-correlation functional by ad¬ 
justment to a dataset of accurate binding energies®®. The 
corresponding parameter is virtually independent of the 
employed reference dataset and hence only depends on 
the employed exchange-correlation functional. 

The independent QHO states can be coupled to the 
long-range part of the Coulomb interaction by solving 
Eqs. (1) and (2) or by explicitly solving the correspond¬ 
ing quasiparticle Hamiltonian®’®®. In the current formu¬ 
lation of MBD the interaction between QHOs is only ac¬ 
counted for via dipole-dipole interactions, enabling an 
analytically exact and efficient solution at vanishing com¬ 
putational expense when compared to DFT. The eigen¬ 
states of this QHO Hamiltonian correspond to collective 
polarization states of the many-body system and the 
corresponding correlation energy is given as the differ¬ 
ence between the zero-point energy of interacting and 
non-interacting QHOs. Therefore, although the initial 
non-interacting response xo is strictly local, the result¬ 
ing interacting response is fully delocalized as would 
be the case in typical ACFDT-based approaches®®’®®. 
In fact, solving the QHO quasiparticle Hamiltonian in 
the dipole approximation is equivalent to calculating the 
ACFDT correlation energy in the RPA for the same set 
of QHOs®®. 

The DFT-bMBD approach has been implemented 
in different software packages including FHI-AIMS®®, 
VASP®^, QuantumEspresso®®, CASTEP®®’®®, ADE®® 
and furthermore is available as an independent code ver¬ 
sion®®. Its success has recently been shown for a va¬ 
riety of systems. Eor small to medium-sized molecules 
the MBD approach correctly captures anisotropy effects 
in the Cq coefficients and the molecular polarizability 
tensor. Eurthermore comparing to experimental dipole- 
oscillator strength distributions this amounts at an ac¬ 
curacy for effective Cq coefficients in DFT-I-MBD of 
about 6.3%, being almost equal to the high accuracy 
achieved in the TS method (5.5%)®. At the same time 
the binding energies are within 5% mean absolute rel¬ 
ative error (MARE) when compared to basis-set limit 
coupled-cluster singles, doubles and perturbative triples 
(CCSD(T)) calculations on the S22 data set of small in- 
termolecular complexes"®. In both cases DFT-I-MBD sig¬ 
nificantly improves on pairwise-additive dispersion ap¬ 
proaches. The strengths of DFT-I-MBD become evi¬ 
dent for extended systems such as large supramolecular 
complexes^®’®® and molecular crystals, where many-body 
interactions play a paramount role. In the latter case, 
pairwise additive approaches fail to achieve the same ac¬ 
curacy they reach for gas-phase intermolecular interac¬ 
tions and overestimate the lattice energy. PBEO-I-MBD 
describes the lattice energy of 16 representative molecu¬ 
lar crystals within a MARE of 4.5% (PBEO-fTS: 12.9% 
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MARE) when compared to lattice energies extrapolated 
from experimental enthalpies of sublimation^^. 

Nevertheless, the dipole approximation utilized in the 
MBD method can sometimes turn out to be insufficient 
and systematic ways exist to extend this approach be¬ 
yond the current state-of-the-art. For example, this can 
be achieved by solving the coupled QHO system with 
an attenuated Coulomb potential. This could be done 
approximately by employing well-established correlation 
techniques from wavefunction theory. In the current work 
we apply the DFT-I-MBD approach for metal-adsorbed 
organic molecules. In the context of dispersion interac¬ 
tions, these are especially challenging systems. On the 
one hand, the localized states of the adsorbate exhibit 
strong attractive van der Waals interactions, on the other 
hand the vanishing band gap of the metal substrate leads 
to a fully non-local collective substrate response that ef¬ 
fectively screens the interactions, thereby reducing Cq 
coefficients, van der Waals radii, and the corresponding 
polarizability^®’®"^. Although this effect can to some ex¬ 
tent be accounted for by effective renormalization of pair¬ 
wise parameters, as has been done in the DFT-l-vdW®'^''^ 
method®®, additional many-body contributions can be 
expected to play a significant role. Although the quan¬ 
tum harmonic oscillators in the MBD calculation have a 
non-vanishing excitation gap and are initially localized, 
the interaction-induced delocalization of the polarizabil¬ 
ity is significantly closer to the correct metallic response 
when compared to a pairwise response. Furthermore, 
the dispersion energy results from an integration over all 
frequencies from 0 to oo [see Eq. (1)]. To supply a good 
starting point for the MBD scheme and to better capture 
the response of the extended substrate, the renormalized 
“atom-in-a-bulk” parameters used in this work are de¬ 
rived using the DFT-|-vdW®““'^ method and subsequently 
employed to parametrize the initial QHO response in the 
MBD scheme. Due to these reasons, we expect the MBD 
method to capture many-body correlation effects with 
reasonable accuracy even in metallic systems. Our ex¬ 
pectation is fully confirmed by the performance of the 
DFT-I-MBD method for realistic molecule-surface sys¬ 
tems discussed in the next section. 


III. RESULTS AND DISCUSSION 

In the following we will study the performance of the 
DFT-I-MBD approach for properties of atoms, molecules, 
and extended nanostructures adsorbed on metal surfaces. 
We study three representative adsorbate-substrate com¬ 
plexes, one of which is dominantly dispersion bound (Xe 
on Ag(lll)), one large molecule adsorbed via both cova¬ 
lent and dispersion interactions (PTCDA on Ag(lll)), 
and an extended organic-inorganic interface (a graphene 
sheet on Ag(lll)). 

All calculations have been performed with the 
DFT-|-vdW®“''^ and DFT-I-MBD implementations in the 
all-electron full-potential FHI-AIMS code using numeri¬ 


cal atomic-orbital basis sets®®. Throughout this work we 
employ the Perdew-Burke-Ernzerhof®® (PBE) and Heyd- 
Scuseria-Ernzerhof®® (HSE) functionals with dimension¬ 
less MBD range-separation parameters of /3 = 0.83 and 
/3 = 0.85 as well as tight numerical basis settings. The 
w range-separation parameter in the HSE functional was 
chosen as 0.11 bohr“^. The binding energy curve for 
Xe on Ag(lll) was performed using the experimentally 
reported (-^3 x •y3)R30° coverage structure for Xe re¬ 
siding at an on-top site of a 6-layered Ag(lll) slab. We 
used a Monkhorst-Pack grid®^ of 15 x 15 x 1 k-points 
in the reciprocal space and a vacuum gap of 20 A. The 
binding energy curve for PTCDA on Ag(lll) was per¬ 
formed using a {- 3 I) surface unit cell in accordance 
to experimental results®®. For reasons of computational 
tractability, the PBE calculations have been performed 
with a 3-layered Ag(lll) metal slab generated with the 
PBE bulk lattice constant, a vacuum gap of 50 A, and a 
Monkhorst-Pack grid of 4 x 4 x 1 k-points in the reciprocal 
space. For the vdW®“''^ and MBD calculations on top of 
PBE, we have used a 5-layered Ag(lll) metal slab to con¬ 
verge the dispersion binding energy. Previous works us¬ 
ing DFT-|-vdW®“‘'^ have shown that geometry relaxations 
induce relatively small deformations for the systems we 
study. In the case of PTCDA on Ag(lll) the change 
in height of the terminal oxygen atoms corresponds to 
0.09 A®®. We therefore assume that this also holds for 
the MBD case. The geometry of graphene has been mod¬ 
eled using a 6-layered slab generated with the PBE bulk 
lattice constant of 4.14 A. It has been optimized using the 
respective dispersion correction method with the bottom¬ 
most layer frozen. In the case of DFT-I-MBD, numer¬ 
ical nuclear forces have been calculated using a finite- 
difference approach. 


A. Xe adsorbed on Ag(lll) 

Noble gas adsorption on metal surfaces has been stud¬ 
ied for a long time and acts as a prototypical example 
for dispersion interactions beyond the simple pairwise 
picture. As shown by Zaremba and Kohn®^, the col¬ 
lective response of the surface modifies the dispersion 
and terms beyond the leading r“® dependence become 
relevant. Taking into account this effect, the vdW®“’'^ 
method projects this interaction into an effective pair¬ 
wise treatment by renormalization of the Cq coefficients, 
atomic polarizabilities, and vdW radii with respect to 
the dielectric response of the substrate, while however, 
still neglecting non-local effects beyond the pairwise ap¬ 
proximation in the treatment of the combined adsorption 
system. 

We have calculated the binding energy curve for Xe on 
Ag(lll) including dispersion interactions with both the 
DFT-|-vdW®“''^ and DFT-I-MBD methods. The adsorp¬ 
tion energy per adsorbed atom was calculated using 

^^ads = A'xe/Ag(lll) ~ (-£'Ag(lll) + Exe) , (4) 
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TABLE I. Adsorption energy, vertical height, and perpendic¬ 
ular vibrational frequency of Xe adsorbed to Ag(lll). Energy 
and vertical height are given for the optimal value as extracted 
from the binding energy curve. For comparison, we also show 
results from experiments. The best estimates are, as given by 
Vidali et 3.6 ± 0.05 A for the adsorption height and an 
adsorption energy between 0.20 and 0.23 eV. 
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FIG. 1. Binding energy curve showing the adsorption energy 
Eads as a function of vertical distance d for Xe on Ag(lll) 
calculated with the HSE-fvdW™’'*' (blue line) and HSE-fMBD 
(green line) methods. The distance d is evaluated with respect 
to the position of the unrelaxed topmost metal layer. The 
results from the binding energy curve and their comparison 
to experiment can be found in Table I. 


where Exe/Ag(iii) is the total energy of the system (Xe 
monolayer -h metal surface), £’Ag(iii) is the energy of the 
bare slab, and Exe is the energy of the isolated Xe gas 
atom, where all quantities correspond to the unrelaxed 
systems. The vertical distance d was defined as the dif¬ 
ference of the position of the atom in the monolayer with 
respect to the position of the unrelaxed topmost metal¬ 
lic layer. Table I shows the optimal adsorption distance 
and energy using PBE and HSE as underlying exchange- 
correlation functionals whereas Figure 1 depicts the bind¬ 
ing curve with the HSE-hvdW"“"f and HSE-hMBD meth¬ 
ods. 

The adsorption height of Xe on Ag(lll) has been stud¬ 
ied experimentally using low-energy electron diffraction 
(LEED)®® and synchrotron x-ray scattering^®, both as¬ 
suming bulk truncation of the metal surface. Adsorp¬ 
tion energetics have been studied using temperature- 
programmed desorption (TPD)^^’^^ and data from in¬ 
elastic helium scattering^®. Including error bars, the ex¬ 
perimentally found adsorption energies and heights range 
from 0.18 to 0.23 eV and 3.45 to 3.68 A. The best esti¬ 
mates, as given by Vidali et are 3.6 ± 0.05 A for 
the adsorption height and an adsorption energy between 
0.20 and 0.23 eV. For comparison, we also show the re¬ 
sults coming from experiments in Table I. 

The results in Table I show an adsorption height of 
3.56 A and an adsorption energy of 0.22 eV with the 
PBE-|-vdW®“'’^ method, demonstrating that already on 
the level of pairwise-additive dispersion including the col¬ 
lective substrate response, the agreement with experi¬ 
ment is excellent. Inclusion of explicit many-body effects 
using the above presented MBD scheme yields, as ex¬ 
pected, only small changes (see Table I). The vertical 


adsorption height is increased above 3.6 A and the ad¬ 
sorption energy is slightly reduced, suggesting a small 
repulsive contribution from higher-order terms beyond 
the leading interatomic behavior. The energy de¬ 
pendence with respect to adsorption height suggests a 
very slow decay with distance from the surface for both 
methods. 

It is clear that a balanced description in terms of ac¬ 
curacy between exchange and correlation can be an es¬ 
sential factor in determining the structural and energetic 
features in adsorption phenomena. Having in mind this 
fact, we have also calculated the binding energy curves 
using HSE as the underlying functional, thereby improv¬ 
ing upon the description of electronic exchange effects. 
These curves are depicted in Figure 1 and the results 
of these are shown in Table I. The adsorption energy is 
not substantially modified when comparing between PBE 
and HSE results, which is not surprising if we consider 
that the attractive part of the interaction in this sys¬ 
tem will mainly be contained in the correlation energy. 
In this sense, the adsorption energy is more sensitive to 
the description of the dispersion interactions than to the 
description of exchange. On the other hand, the adsorp¬ 
tion height seems to be more sensitive to the choice of 
the exchange-correlation functional. The vertical adsorp¬ 
tion height is 3.52 A with the HSE-|-vdW®“"'^ method and 
3.57 A with the HSE-I-MBD method. Inclusion of many- 
body effects using the MBD scheme yields a slight in¬ 
crement of 0.05 A in the adsorption height when HSE is 
the underlying exchange-correlation functional, in con¬ 
trast to the larger increment of 0.08 A observed with 
the PBE functional. These facts suggest that the inclu¬ 
sion of screened short-range exact exchange, as found in 
the HSE functional, provides a more balanced descrip¬ 
tion when coupled to the MBD scheme as part of the 
correlation energy. Although the effects may seem to be 
small in the case of Xe on Ag(lll), these may become 
increasingly important in the adsorption of organic ad¬ 
sorbates on metal surfaces, where we find a much more 
complex interplay of interactions. As a final remark, it is 
worth to mention that many-body effects persist for large 
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distances (larger than approximately 5.0 A) as the com¬ 
parison between binding curves with the vdW’*“''^ and the 
MBD methods in Figure 1 shows. The discrepancy be¬ 
tween the binding curves at an infinite distance is given, 
in accordance to the definition in Eq. (4), by the dif¬ 
ferent values for the formation energy of the monolayer 
calculated with the vdW®”^ and the MBD methods. 

Rohlfing and Bredow have studied Xe on Ag(lll) us¬ 
ing a correlation treatment based on the RPA including 
exact exchange (EXX)*^, which accounts for many-body 
effects by explicitly calculating the long-range correla¬ 
tion energy of the system. The corresponding adsorp¬ 
tion energy lies about 30 meV above our PBE-I-MBD 
and HSE-I-MBD results, whereas the adsorption height 
is found to be in good agreement with both experiment 
and only 0.04 A below PBE+MBD and 0.03 A above our 
HSE+MBD results. This suggests that explicit account 
of many-body dispersion does in fact reduce the bind¬ 
ing strength in comparison to a pairwise treatment. In 
the latter case of correlation based on the RPA in Ref.^, 
this effect might be overestimated due to neglect of the 
exchange-correlation kernel, the underlying plasmon-pole 
approximation, and the fact that the response function 
of the system is not fully coupled, being calculated sepa¬ 
rately for substrate and adsorbate. 

We have also computed the perpendicular vibrational 
frequencies of Xe on the metal surfaces to probe the 
curvature of the potential energy curves around the 
minimum in each case. Eor this, following previous 
works^®^^®, we have modeled the gas-surface adsorption 
potential with the following function given by the sum of 
repulsive and attractive dispersion interactions 

EMd) = + Aml, (5) 

[d- ZoY 

where Aads (d) is the adsorption potential between Xe and 
the metal substrate at a distance d from the surface and 
Eml is a constant that corresponds approximately to the 
formation energy of the Xe monolayer. We have deter¬ 
mined the parameters oi, a 2 , Ca, Zq, and Eml by fit¬ 
ting Eq. (5) to the binding energy curve calculated with 
each method. The vibrational energy Ayib is then cal¬ 
culated using Evih = hv = hj{2 tt) where v, 

h, and mxe are the vibrational frequency, Planck’s con¬ 
stant, and the mass of a Xe atom, respectively. The force 
constant ke corresponds to the second derivative evalu¬ 
ated at the minimum of the potential given by Eq. (5). 
Eollowing this procedure, the results for Avib are given 
in Table I. Considering the absolute magnitude of this 
vibration known from experiment (2.8 meV) all methods 
perform very well, with explicit many-body dispersion re¬ 
ducing the vibrational energy closer towards the experi¬ 
mental value in the cases of PBE-I-MBD and HSE-I-MBD. 
This MBD-induced reduction in frequency is also visible 
from the width and curvature of the HSE-I-MBD binding 
energy in Figure 1 when compared to HSE-|-vdW®“‘'^. 

Although we are naturally interested in the computa¬ 
tion of adsorption distances and energies and how they 



FIG. 2. Differential static polarizability Aa® of Xe on Ag(lll) 
[see Eq. (6)] as a function of the vertical distance d of the Xe 
monolayer as calculated with the HSE-I-MBD method. Shown 
are the components parallel (blue curve) and perpendicular 
(red curve) to the surface plane. 

compare to experiments, we must also understand the 
improvements in the description of the physics and chem¬ 
istry behind our methods. One of the novelties behind 
the MBD method is a more robust description for the po¬ 
larizability of molecules and solids as the method includes 
many-body screening effects coming from the electrody¬ 
namic response of the system, a fact which leads to a 
more accurate description of dispersion interactions (see 
section H). A more detailed understanding of these ef¬ 
fects in the adsorption properties of Xe on Ag(lll) can 
be gained by studying the changes in the static polar¬ 
izability tensor of the system Aq:° as a function of the 
adsorption height d of the Xe monolayer. We have com¬ 
puted these changes using 

Ad^ll/_L) = d(||/±)gy^ - («(llA)Ag(lll) + “°iA)xe) ’ 

( 6 ) 

where is the screened static polarizability of 

the complete system (Xe monolayer -|- metal surface), 
®Ag(iii) screened static polarizability of the bare 

slab, and screened static polarizability of the 

Xe monolayer in periodic boundary conditions. The 
treatment of the dipole-dipole coupling between atoms in 
the MBD method naturally introduces anisotropy in the 
polarizability of the complete sytem®. Due to the symme¬ 
try of the system, we observe two equivalent components 
in the polarizability of the system which lie parallel to the 
plane of the surface; we denote their change as Ad®. The 
third component points in the direction perpendicular to 
the plane of the surface, we denote its change as Ad° . 
Figure 2 shows the results of Eq. (6) for each component 
as a function of the adsorption height d. The quantities 
shown in Eigure 2 correspond to the MBD calculations 
associated to each distance in the binding energy curve 
of Figure 1. 
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FIG. 3. PTCDA overlayer adsorbed to Ag(lll) surface with 
the supercell geometry shown as a dotted line. 


The definition of Aq:° in Eq. (6) lets us identify the 
distance upon adsorption at which the coupling of the 
components becomes relevant for the polarizability of the 
system. At large adsorption distances, the polarizability 
of the complete system, given by the first term in Eq. 
(6), is equal to the sum of its parts, found in the second 
term (in parenthesis) of Eq. (6). This balance towards 
zero starts approximately at 5.5 A and becomes practi¬ 
cally zero at distances greater than 7.5 A no matter which 
direction is considered. At distances lower than 4.5 A, 
the coupling between Xe and Ag(lll) starts to become 
relevant. Upon reduction of the adsorption distance of 
the Xe monolayer, the polarization of the system in the 
direction parallel to the surface (given by the blue curve) 
is decreased in favor of an increasing polarization in the 
direction perpendicular to the surface (given by the red 
curve) due to the interaction with the surface. This is re¬ 
flected in the negative values found in Ady upon adsorp¬ 
tion as well as the positive values in Ad*|. This behavior 
yields a stronger polarization of the system in the direc¬ 
tion perpendicular to the surface plane at the equilibrium 
adsorption distance. Correspondingly, the increased po¬ 
larizability towards the surface leads to larger interac¬ 
tion screening from many-body contributions, while at 
the same time dispersion interactions between Xe atoms 
are reduced. The changes in the anisotropic terms of 
the polarizability of the system may seem to be small 
in magnitude. However, these subtle changes induced 
by the coupling of the complete system could generate a 
preferential interaction along a specific direction parallel 
or perpendicular to the surface plane, yielding direction¬ 
ality in the formation of molecular monolayers. 


B. PTCDA adsorbed on Ag(lll) 


3,4,9,10-perylene-tetracarboxylic acid (PTCDA) is 
one of the best-studied large organic adsorbates on 
coinage metal surfaces, both experimentally^®^®^ and 
theoretically^’^^’®^’®^. PTCDA consists of two anhydride 
head-groups and an aromatic backbone. The molecule 


adsorbs in a densely-packed herringbone arrangement on 
the Ag(lll) surface (see Figure 3), with the oxygens be¬ 
ing tightly bound via covalent interactions and the aro¬ 
matic backbone interacting with the surface via disper¬ 
sion interactions^®. Using current electronic structure 
methods, either based on pairwise dispersion-corrections 
or RPA-based approaches, the interaction energy as well 
as the vertical adsorption height of PTCDA adsorbed on 
the Ag(lll) surface is on the upper range of what is ex¬ 
pected to be the experimental adsorption energy (vide 
infra). In the case of PBE-l-vdW®'^'’^, this happens even 
though the metallic substrate response has already been 
effectively accounted for in the dispersion parameters. 
As it has also been shown for the case of azobenzene 
adsorbed to Ag(lll), neglecting this effect leads to even 
further overestimation of adsorption energies and adsorp¬ 
tion heights®^. The remaining deviations to experiment, 
for example the overestimation of the adsorption interac¬ 
tion, may be ascribed to missing many-body dispersion 
contributions and the semi-local treatment of exchange 
interactions, both of which we will discuss below. 

On the experimental side, the adsorption geometry and 
the individual atomic adsorption heights of PTCDA on 
Ag(lll) are known from x-ray standing wave (XSW) 
measurements®®. The adsorption energy, typically ex¬ 
tracted from TPD measurements, however, is not ex¬ 
perimentally known due to the molecule being destroyed 
upon thermal desorption®^. However, Stahl et al. have 
measured a binding energy of 1.16 ± 0.1 eV for 1,4,5,8- 
naphthalene-tetracarboxylic-dianhydride (NTCDA) on 
Ag(lll)®®, a smaller analogue of PTCDA containing 
the same number of terminal oxygen atoms but with 
a smaller aromatic backbone (40% less molecule sur¬ 
face area). Based on surface-state photoemission data, 
Ziroff et al. estimated the binding energy of PTCDA 
and its smaller analogue NTCDA on Au(lll) by as¬ 
suming purely physisorptive surface-binding and utiliz¬ 
ing a connection to noble-gas adsorption on metals. The 
corresponding binding energies are reported as 2.0 and 
1.5 eV for PTCDA and NTCDA on Au(lll)®®. Con¬ 
trary to that, Wagner et al. report significantly higher 
binding energies of 2.6 and 1.7 eV for both molecules on 
Au(lll) from force pulling experiments using an atomic 
force microscope (AFM)®^. Considering the TPD data 
on NTCDA on Ag(lll) and the fact that two experi¬ 
ments place the additional binding energy of PTCDA on 
Au(lll) compared to NTCDA on Au(lll) to be within 
33 to 66%, we estimate the adsorption energy of PTCDA 
on Ag(lll) to be between 1.4 and 2.1 eV (see Table H). 

We have calculated the binding energy curve for a 
PTCDA monolayer on Ag(lll) (two adsorbed molecules 
per unit cell, see Figure 3) including dispersion inter¬ 
actions with both the PBE-l-vdW®'^’'^ and PBE-I-MBD 
methods. The adsorption energy per molecule in the 
monolayer was calculated using 

Eads = 2 [^PTCDA/Ag(lll) ~ (-E'Ag(lll) +-UptCDa)] , 

( 7 ) 













TABLE II. Adsorption energy and vertical height of PTCDA 
adsorbed to Ag(lll). Energy and vertical height are given for 
the optimal value as extracted from the binding energy curve. 
For comparison, we also show results from experiments. 


PTCDA/Ag(lll) 

-Uada [eV] 

d[M 

PBE-tvOW^^"*" 

2.50 

2.89 

PBE-tMBD 

1.77 

2.94 

vdW-DF“°“-"‘=®3 

2.0 

3.5 

vdW-DF-cx^'^ 

3.5 

3.1 

cRPA-tEXX'^ 

2.4 

3.1 

Exp.80.85-87 

(1.4- 2.1)“ 

2.86 ± 0.01 


“ estimate from experimental results: see text for more 
details 
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FIG. 4. Binding energy curve showing the adsorption en¬ 
ergy i?ads as a function of vertical distance d for PTCDA on 
Ag(lll) calculated with the PBE-bvdW"”*' (blue line) and 
PBE-I-MBD (green line) methods. The distance d is evalu¬ 
ated with respect to the position of the unrelaxed topmost 
metal layer. The results from the binding energy curve and 
their comparison to experiments can be found in Table II. 

where i?pTCDA/Ag(iii) is the total energy of the system 
(PTCDA monolayer -b metal surface), i?Ag(iii) is the 
energy of the bare slab, and Aptcda is the energy of 
the PTCDA monolayer in periodic boundary conditions, 
where all quantities correspond to the unrelaxed systems. 
The vertical distance d was defined as the difference of 
the position of the monolayer with respect to the unre¬ 
laxed topmost metallic layer. Table II shows the optimal 
adsorption distance and energy using PBE as underlying 
exchange-correlation functional whereas Figure 4 depicts 
the binding energy curve with both methods. Includ¬ 
ing many-body dispersion contributions via PBE-I-MBD, 
we find that the adsorption binding strength is reduced 
in comparison to PBE-bvdW®'^'^^. This is reflected in a 
reduced interaction energy and an increased adsorption 
height as depicted in Figure 4. Accounting for the higher- 
order correlation terms and the correct intermolecular 


polarization counteracts the dispersion energy of individ¬ 
ual pairs of atoms stemming from the leading term. 
From the binding-energy curves in Figure 4 we further 
find that the curvature around the basin of attraction is 
reduced with the basin itself being widened. As in the 
case of Xe on Ag(lll) this suggests both a reduced vi¬ 
brational frequency and larger anharmonic contributions 
along the molecule-surface mode. The higher-order dis¬ 
persion terms from the MBD come into effect at practi¬ 
cally all distances considered in Figure 4 reducing the dis¬ 
persion energy. The MBD binding energy closely follows 
the vdW®'^“'^ binding energy only at adsorption distances 
larger than approximately 6.0 A. Upon further reduction 
of the adsorbate-surface distance, the electron density 
overlap is increased and the Pauli repulsion becomes the 
dominant term. At the same time, the many-body contri¬ 
butions become smaller by approaching a distance from 
the first metal layer of the order of the range-separation 
length. 

By inclusion of explicit many-body dispersion we also 
find an improvement in performance when compared to 
experiment and literature data. Whereas PBE-l-vdW®'^'’^ 
is known to already yield a good description of ad¬ 
sorbate geometries^®’^^, corresponding adsorption en¬ 
ergies are systematically overestimated (see Table II). 
The PBE-I-MBD scheme yields an adsorption energy of 
1.77 eV that lies within the estimated regime, and an 
adsorption height of 2.94 A which exceeds the experi¬ 
mental adsorption distance of 2.86 A by 0.08 A. Com¬ 
paring to other calculations based on explicit correlation 
treatment, either on the RPA' or the vdW-DF leveU"^’®^, 
we find that PBE-I-MBD yields lower adsorption energies 
at smaller equilibrium distances from the surface. The 
overestimation of equilibrium adsorption distances is a 
well known issue for first and second generation vdW- 
functionals such as vdW-DF^^’^^. However, this prob¬ 
lem seems to have been remedied to some extent in the 
most recent vdW-DF-cx functionaU^. In the case of the 
PBE-I-MBD method, the adsorption distance is within 
0.1 A of the value found in experiment. Moreover, an 
accurate treatment of screened exchange in the under¬ 
lying functional, via e.g. HSE, will further improve the 
description of the geometry as happens in the case of Xe 
on Ag(lll), where we have observed a reduction of the 
adsorption distance by 0.07 A upon coupling the MBD 
energy to the HSE functional (see section HI A). Finally, 
it is important to point out that our calculations for the 
binding-energy curve correspond to an unrelaxed system, 
that is a planar PTCDA monolayer adsorbed on an un¬ 
relaxed surface slab. Undertaking a full relaxation of 
the system using the MBD method will lead to the well 
known distortion of the molecule within the monolayer 
thereby yielding an improved description of the adsorp¬ 
tion distance. For example, full geometry relaxation of 
PTCDA/Ag(lll) using the PBE-|-vdW®“''^ method de¬ 
creases the average adsorption distance by 0.09 A and 
increases the binding energy by 0.36 eV. The same change 
with PBE-I-MBD would still place the results within the 
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FIG. 5. Graphene sheet adsorbed to Ag(lll) in a \/Z x \/3 
surface unit cell. 


expected experimental range. 


C. Graphene adsorbed on Ag(lll) 

Graphene monolayers adsorbed on metal surfaces rep¬ 
resent especially challenging systems with regards to dis¬ 
persion interactions. Graphene sheets and more generally 
carbon nanostructures, such as carbon nanotubes and 
fullerenes, exhibit different dispersion interactions than 
small organic molecules. This is apparent from the signif¬ 
icant deviations of effective Cg coefficients across different 
carbon-based materials when accounting for polarization- 
induced screening®®. The reason for this is a largely delo¬ 
calized, collective polarization response to electric fields. 
Equally for a graphene monolayer interacting with an 
Ag(lll) surface, screening effects and many-body dis¬ 
persion contributions to the adsorption energy can be 
expected to be large. We have performed calculations for 
a graphene sheet commensurately adsorbed on a y/S x ^/S 
Ag(lll) slab (see Fig. 5). In this adsorption geom¬ 
etry, the intramolecular carbon bonds in graphene are 
slightly elongated with 2.54 A, when compared to the 
2.46 A of isolated graphene, however the unit cell is still 
small enough not to induce significant buckling in the 
graphene sheet due to the metal surface corrugation. We 
choose this unit cell for reasons of comparability to other 
studies^^’®® and computational tractability. We define 
the adsorption energy of graphene (Gr) per carbon atom 
as 

Eads = [EGr/Ag(lll) ~ (EAg(lll) + Ecr)] /Ec. (8) 

The PBE-I-MBD approach reduces the adsorption en¬ 
ergy per atom of graphene on Ag(lll) over 38% com¬ 
pared to the pairwise PBE-|-vdW®““'^ approach. The re¬ 
sulting 45 meV per carbon atom are found at an equilib¬ 
rium adsorption height of 3.23 A. As found for Xe and 
PTCDA on Ag(lll), the equilibrium adsorption height 
is increased by inclusion of many-body effects. However 


the effect on the adsorption height can be considered 
strong when compared to Xe and PTCDA and stems 
from the large magnitude of many-body effects between 
metal and graphene that effectively screen the dispersion 
interactions. The resulting binding energy of 45 meV 
per carbon atom is close to the interlayer binding en¬ 
ergy of crystalline graphite as predicted by PBE-I-MBD 
(48 meV/C atom)®®, suggesting a consistent description 
of many-body effects for nanostructured carbon. 

The importance of many-body contributions is evident 
from the collective eigenvectors of the QHO Hamiltonian 
that represent the eigenstates of the long-range correla¬ 
tion problem in the basis of atomic positions. Although 
these eigenstates merely constitute the canonical basis of 
the QHO model Hamiltonian and do not correspond to 
any actual physically observable quantities, their change 
upon adsorption can help to qualitatively analyze the col¬ 
lective mutual polarization and depolarization between 
different domains of the system. Figure 6 shows two 
such representative MBD eigenmodes (top and bottom) 
for an isolated graphene sheet and graphene adsorbed on 
Ag(lll) (left and right). The blue vectors indicate the 
direction and magnitude of polarization on each localized 
QHO (depicted as orange spheres). The first eigenmode 
(a and b in Fig. 6) describes lateral polarization within 
the graphene sheet and is almost unaffected upon ad¬ 
sorption on the metal surface. The second mode (c and 
d in Fig. 6) represents polarization orthogonal to the 
graphene plane and is strongly modified due to adsorp¬ 
tion. In fact, this mode is fully delocalized over adsorbate 
and substrate (not shown here) and describes the collec¬ 
tive polarization between the subsystems. These visually 
apparent changes can also be seen in the energy of the 
eigenstates. Whereas the energy of the fist eigenmode 
only decreases by about 1.7 meV upon adsorption, the 
second mode contributes about 32 meV to the dispersion 
energy defined by the sum of eigenenergies. Due to the 
energy shifts of the MBD modes we can in fact pinpoint 
which polarization modes yield the most important con¬ 
tributions to the dispersion energy. 

It is possible to visualize the change in MBD modes 
due to adsorption analogous to a density-of-states (DOS) 
diagram (see Fig. 7). In the MBD-DOS of graphene ad¬ 
sorbed to Ag(lll), the first leading peak almost solely 
stems from polarization modes of the graphene sheet as 
apparent from the comparison to the DOS of isolated 
graphene (blue curve in Fig. 7), whereas the second large 
peak almost solely stems from eigenenergies of the sub¬ 
strate. The actual modifications in the graphene DOS 
due to adsorption can be seen by projecting the iso¬ 
lated graphene states from the spectrum of graphene 
on Ag(lll) (red curve). Whereas the overall DOS of 
graphene appears basically unchanged, a number of new 
modes appears at higher energies around 16 eV due to 
hybridization with modes of the substrate. At the same 
time polarization modes previously localized to the sub¬ 
strate are shifted towards lower energies (magenta curve) 
by hybdridization with graphene (at ~ 11.5 eV in Fig. 
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FIG. 6. Visualization of MBD eigenmodes for graphene on Ag(lll). Orange spheres represent QHOs at the positions of the 
atoms, arrows depict the polarization direction. Shown are a lateral polarization mode within the graphene sheet for the isolated 
(a) and adsorbed sheet (b) and a vertical polarization mode, again, for the isolated (c) and adsorbed sheet (d). Contributions 
from the underlying metal surface are omitted. 



MBD eigenstate energy in (eV) MBD eigenstate energy in (eV) MBD eigenstate energy in (eV) 


FIG. 7. Quasiparticle Density-of-State (DOS) of MBD eigenstates for graphene adsorbed on Ag(lll) (left), isolated graphene 
and the projected DOS of graphene (center, red and blue curves), and a clean Ag(lll) surface and the projected Ag(lll) DOS 
(right, green and magenta curves) contributions. 


7). While these changes in eigenmodes are instructive to 
study, it should be noted that they correspond to fluc¬ 
tuations of a model system composed of coupled QHOs. 
Therefore, unfortunately comparison of the MBD modes 
to experimental observables is not possible. Only inte¬ 
grated quantities in the MBD model, such as frequency- 
dependent polarizability or the binding energy, can be 
directly compared to experiment. 

The MBD framework enables us to capture the qual¬ 
itative physics of many-body dispersion and interpret 
it with concepts familiar to electronic structure the¬ 
ory. On the other hand the quantitative performance of 
PBE-I-MBD, in the case of graphene on Ag(lll), can only 


be evaluated in comparison to experimental data. How¬ 
ever, due to a lack of such data for graphene on Ag(lll), 
comparison to other simulation approaches can help to 
put binding energies and adsorption geometries into per¬ 
spective. While the effective pairwise PBE-|-vdW®“’'^ 
approach already improves on the notorious underbind¬ 
ing of pure semi-local functionals, the resulting adsorp¬ 
tion energy and height are still above and below what 
is found when applying non-local van der Waals func- 
tionals®®^®^ or RPA-based correlation methods^^ (see Ta¬ 
ble IB). Introducing many-body effects via MBD puts 
the resulting adsorption height and adsorption energy 
in the same range as for example optimized van der 
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TABLE III. Adsorption energy and vertical height of a 
graphene sheet adsorbed to Ag(lll). Energy is given in meV 
per carbon atom, the distance is measured as average distance 
between carbon atoms and the first substrate layer. 


Graphene / Ag( 111) 

-Aads/C atom 

[meV] d [A] 

PBE-tvdW^^"' 

72 

3.05 

PBE-tMBD 

45 

3.23 

LDA**^*’ 

30 

3.22 

RPBE*^° 

1 

5.57 

(cRPA-tEXX)@PBE2° 

78 

3.31 

vdW-DF2(C09x)«® 

58 

3.27 

optB86b-vdW-DF®® 

67 

3.34 

vdW-DF®** 

33 

3.55 

rVVlO®® 

68 

3.48 


Waals functionals as proposed by Klimes and Michaelides 
(optB86b-vdW-DF)^^’®®. Results obtained from exact 
exchange and correlation based on RPA (cRPA+EXX)^° 
yield an adsorption height similar to what is found with 
PBE+MBD, but at a surprisingly high adsorption energy 
per atom, considering that cRPA is typically considered 
to underestimate binding energies of small intermolecu- 
lar complexes®^. In fact, the 78 meV/ C atom adsorption 
energy reported by Olsen et al. might not be fully con¬ 
verged with respect to the number of substrate layers, 
basis-set truncation, and k-point sampling, considering 
the slow convergence of RPA correlation energies with 
respect to these parameters'*^. Additionally, the neglect 
of the q = 0 contribution to the correlation energy might 
yield an unexpected bias towards higher interaction en¬ 
ergies. The overall remaining spread in the adsorption 
energy and height for graphene/Ag(lll) across methods 
still poses a challenge to be further investigated, for which 
experimental reference values would be welcome. 


IV. CONCLUSIONS 

Accurately treating electronic correlations in realistic 
systems without resorting to highly demanding corre¬ 
lation methods, such as the ’chemistry gold standard’ 
CCSD(T)®'^ or ACFDT-based correlation techniques, can 
be exceedingly challenging. This is even more so the case 
for condensed-matter systems, where periodicity, simul¬ 
taneous existence of localized and delocalized states, and 
the shear system size, limit the applicability of quantum- 
chemical approaches. The here discussed DFT-I-MBD 
method yields an explicit account of long-range many- 
body dispersion effects, including instantaneous polar¬ 
ization effects, collective response or screening, which is 
not only necessary to yield a quantitative account of ad¬ 
sorption energies, but often also is necessary to yield a 
qualitatively correct description of adsorbate-substrate 
structure and binding. In this work, we analyzed the 
importance of many-body dispersion effects for the cor¬ 


rect description of adsorption on metal surfaces for the 
three test cases of Xe, PTCDA, and graphene adsorbed 
on Ag(lll). We furthermore evaluated the accuracy and 
performance of DFT-I-MBD in capturing these effects 
compared to other approaches and experiment. 

Although many-body contributions did not strongly al¬ 
ter the interaction strength for noble-gas atom adsorption 
in the case of Xe on Ag(lll), the effects were still signifi¬ 
cant considering the system size and furthermore consid¬ 
ering the role of noble-gas atoms in the development and 
justification of pairwise dispersion-correction approaches. 
Most interesting is the strong anisotropy of the atomic 
polarizability tensor upon adsorption that is captured 
with MBD. In the cases of PTCDA and graphene, inclu¬ 
sion of many-body effects reduces the overbinding of pair¬ 
wise dispersion by 20 and 27 meV per adsorbate atom, 
amounting to a reduction of adsorption energy of approx¬ 
imately 29% and 25% for both systems. At the same time 
the adsorption height was increased by about 0.1 and 
0.2 A, respectively. This shows that a mere account of 
pairwise additive dispersion leads to significant overbind¬ 
ing and many-body effects must not be neglected. The ef¬ 
fect of higher-order many-body dispersion contributions 
is directly visible in the MBD eigenstates that describe 
the delocalized polarization between adsorbate and sub¬ 
strate. The subsystem hybridization visible from the 
eigenmode analysis can function as an interpretational 
tool to understand the adsorbate-substrate interactions. 

Despite the clear success of the DFT-I-MBD approach, 
several methodological (c/. Ref. 5) as well as issues spe¬ 
cific to metal-surface adsorption remain. One of which 
is that the polarizability screening and many-body in¬ 
teractions described by the DFT-I-MBD method do not 
yet explicitly feed back into the electronic density or 
affect the potential that enters the Kohn-Sham equa¬ 
tions. As a result the DFT level alignment of adsor¬ 
bate states with respect to the substrate Fermi level will 
not be directly affected by the MBD correlation con¬ 
tribution and will still be described at the level of the 
semi-local exchange-correlation treatment. This effect 
together with the self-interaction error in the exchange 
functional will still plague short-range interactions such 
as covalent bonding. However, in the above presented 
range-separation framework, density functionals beyond 
the generalized gradient approximation can also be cou¬ 
pled with the MBD scheme, as was already shown above 
with the application of HSE-I-MBD for Xe and PTCDA 
on Ag(lll). 

As was shown in section II of this manuscript, in treat¬ 
ing adsorbates on metal surfaces DFT-I-MBD can be un¬ 
derstood as combining, so to say, the best of both worlds 
- of DFT and wavefunction approaches. Whereas the 
semi-local or screened hybrid density functionals cor¬ 
rectly treat the delocalized metallic states of the under¬ 
lying substrate, the missing long-range correlation is ac¬ 
counted for by recasting the long-range correlation prob¬ 
lem into an auxiliary system of coupled QHOs. The close 
ties of the MBD approach to correlation techniques open 
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a path to systematically go beyond the current state- 
of-the-art. Replacing the dipole-dipole interactions with 
an attenuated Coulomb potential, the full machinery of 
quantum chemistry could be employed to solve the cou¬ 
pled QHO Hamiltonian. The currently employed “min¬ 
imal basis” approach of associating every atom with a 
single QHO could furthermore be extended by an expan¬ 
sion of the electron density around an atom in terms of a 
linear combination of QHOs. With such a basis set, the 
coupled QHO problem can be solved using any correlated 
quantum-chemical method, all of which would still be at 
a comparably low computational cost when compared to 
solving the full electronic many-body problem. 

The above results show that the DFT-I-MBD approach, 
although solely developed for describing the correlation 
effects in molecules and finite-band gap materials, im¬ 
proves the description of molecular adsorption on metal 
surfaces for the here studied systems. The initially local¬ 
ized QHOs do not correctly account for the response of 
the metal, however, their interaction leads to a delocal¬ 
ized response that significantly improves the description 
of long-range correlation. As a result the description of 
the density-density response becomes closer to a metallic 
response. In addition to a gain in accuracy, including 
many-body dispersion effects via the MBD technique of¬ 
fers many conceptual insights into adsorbate-substrate 
bonding, which, overall, makes it a strong contender in 
the accurate modeling of complex organic-inorganic in¬ 
terfaces. 
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